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Abstract 

We investigate dynamics of a self-propelled de- 
formable particle under external field in two dimen- 
sions based on the model equations for the center of 
mass and a tensor variable characterizing deforma- 
tions. We consider two kinds of external force. One 
is a gravitational-like force which enters additively in 
the time-evolution equation for the center of mass. 
The other is an electric-like force supposing that a 
dipole moment is induced in the particle. This force 
is added to the equation for the deformation tensor. 
It is shown that a rich variety of dynamics appears by 
changing the strength of the forces and the migration 
velocity of self-propelled particle. 



1 Introduction 

Dynamics of self-propelled objects have attracted 
much attention recently from the view point of var- 
ious new orders far from equilibrium. In order to 
elucidate the new dynamical states, one needs to 
develop the theory of nonlinear dynamics and non- 
equilibrium statistical physics. There are many kinds 
of self-propelled motions both in biological and non- 
biological systems. Migration of microorganisms such 
as living cells and bacteria [TJ [2] [3] [4] [5] [6] are typi- 
cal examples in biological systems. Artificial models 
have been introduced for self-organized microswim- 
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mers [7] [8] [9]. On the other hand, there are investiga- 
tions of self-propulsion in non-biological systems such 
as oily droplets in surfactant solution [TQl [11] , colloid 
composites under external field or with chemical re- 
actions [HI [131 [14] and Janus particles pT5l [TBI [17] . 
Marangoni effects due to chemical reactions and elec- 
trophoresis are the origins of the self-propulsion of 
the above systems. See a review article [18] for the 
related nano-systems. It is known that theoretical 
formulation of self-propulsion is non-trivial even for 
a single particle. 

In a previous paper, one of the present authors 
(TO) and Ohkuma introduced a set of model equa- 
tions for a deformable self-propelled particle [19] , 
The model equations have been constructed by sym- 
metry argument and hence from a general view point 
of a dynamical system. After that, those equations 
together with higher modes of deformations have 
been derived by the singular perturbation method 
for an excitable reaction-diffusion system [2UJ both 
in two dimensions [3T] and in three dimensions [22] . 
Various dynamics of motions are obtained by chang- 
ing the migration velocity and the softness of the 
particle [19] [23] [24] . The collective dynamics of the 
particles with oricntational interactions have been in- 
vestigated in two dimensions [T21 US US] . It is men- 
tioned that deformable self-propulsion was studied by 
Shapere and Wiczek to investigate efficiency of swim- 
ming motion |27] . It is also mentioned that there are 
a number of theoretical studies of collective dynam- 
ics of non-deformable self-propelled particles. Sec 
Ref. [2S] and the earlier papers cited therein. 

In the present paper, we study dynamics of a de- 
formable self-propelled particle under external force 
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in two dimensions. The system we consider is a sin- 
gle isolated particle but has internal degrees of free- 
dom due to deformability. It should be noted that, 
even when the external forcing is absent, there is a 
bifurcation between a straight motion and a circu- 
lar motion along a closed trajectory [T5]. Therefore, 
by adding an external force, there occurs a conflict 
or frustration between the circular motion and the 
forced straight motion. As a result, non-trivial be- 
havior of motion and bifurcations are exhibited by 
changing the magnitude of the external force. A 
model of self-propulsion which undergoes a circular 
motion has also been introduced in Rcf . 29 . 

We consider two kinds of external forcing. One is a 
gravitational-like force which enter additively in the 
equation of motion for the center of mass. A charged 
particle under electric field is also governed by this 
euation. The other is an electric-like force which is 
supposed to produce an electric dipole in the particle. 
The force is added to the equation for the deforma- 
tion tensor. We shall carry out numerical simulations 
of the set of the time-evolution equations to obtain a 
dynamical phase diagram. Analytical study will also 
be developed to reproduce some of the motions and 
the bifurcations. 

The organization of this paper is as follows. In 
the next section (section [2]) , we introduce the model 
equation. Numerical results for the gravitational-like 
forcing is given in section [3J We will show that a 
circular-drift motion occurs in this case, where a par- 
ticle causes a drift motion to the direction perpendic- 
ular to the external force. This motion is analyzed 
in detail in section [4] Numerical simulations for the 
electric-like forcing are described in section [5] An- 
alytical study of the gravitational-like force is pre- 
sented in section [6] whereas that of the electric-like 
force is given in section Summary and discussion 
are given in section [8] 



2 Model equations 

We consider the following set of equations of motion 
in two dimensions for the velocity of the center of 
mass v = (yi, v^) and the symmetric tensor S a p char- 



acterising deformations 



dv, 
~dt 



j- = iv a - \v 2 \v a - aS a0 vp + g a (1) 
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This type of external force has been considered in a 
phase separated droplet under electric field [30] . The 
coefficient 7 may change the sign and k is positive. 
a and b are coupling constants. The tensor S a p is 
traceless to make the area of the particle (volume of 
the particle in three dimensions) constant and defined 
by S a p — s (n a np — Tj6 a ,p), where the unit vector 
n is parallel to the long axis of a deformed elliptical 
particle and s > is the degree of deformation from a 
circular shape. The external force g a is assumed to be 
given by g = (0, —g) with g > 0. The other external 
force E a in Q a p is applied as E = (1,0). The left- 
right and the up-down symmetry of the system allow 
us to put h > for fixed signs of a and b. The set 
of equations (p} and ([2]) without the external forces 
has been introduced in Ref. [19] and has been derived 
from an excitable reaction-diffusion system not only 
in two dimensions but also in three dimensions |21[ 
[22]. 

Equations (jTJ and ([2]) can be written as 
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with v and s positive values, and 



Putting Eqs. (0 and (O together, we have 



(12) 



dt 



1 fbv 2 

— as sin 2t/> H — cos < 

2 \ s / u 



2s 



sin2(0 + V). 



(13) 



Equations ([T]) and ([2]) in the absence of g and 
Qap — exhibit a drift bifurcation. That is, when 
7 < 0, a stable state is a motionless state, and when 
7 > 0, a particle undergoes a self-propelled motion. 
In the latter case, there is another bifurcation [T9"] . 
The constant defined by 



2k 



(14) 



plays an important role. If B is positive as we assume 
throughout this paper, a particle undergoes a straight 
motion in some direction determined by the initial 
condition for < 7 < j c . The bifurcation threshold 
7 C is given by [19] 



7c 



K 

ab 



K _ k(1 + B) 
2 ~~ 



2B 
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There are two cases of the straight motion depending 
on the sign of b. If b is positive, the particle elongates 
along the direction of the migration velocity whereas, 
if it is negative, the elongation is perpendicular to the 
velocity. This straight motion loses its stability for 
7 > 7 C and a periodic motion along a closed circle 
appears to be stable |19) . 

In the following sections, we study the interplay 
between the circular motion and the straight motion 
forced by the external fields. 

3 Numerical Results I 

In this section, we show the results of numerical 
simulations of the model equations ((4]) - (JT]) for the 
gravitational-like external force, i.e., with Q a p = 0. 
The fourth-order Runge-Kutta method is employed 
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Figure 1: (Colour on-line) Phase diagram on the 7 — <7 
plane for (a) k = 0.2 and (b) k = 0.75. The sym- 
bols indicate the following motions; the circular-drift 
motion (cross), the zigzag- 1 motion (diamond), the 
zigzag-2 motion (square), and the straight-falling mo- 
tion (circle). Note that there are coexistence regions. 
The thin solid line is the saddle homoclinic-orbit bi- 
furcation boundary determined numerically from the 
reduced equations (| and (Ti~9"|) . The thick solid 
line is the Hopf bifurcation boundary obtained by 
Eq. (j45|). 



with time increment St = 10~ 4 . The coupling coef- 
ficients a and b are fixed as a = —1.0 and b = —0.5. 
We have obtained the phase diagram on the j-g plane 
for k — 0.2 and K = 0.75, as shown in Fig. [T] The bi- 
furcation threshold is given by 7 C = 0.18 for n = 0.2 
and 7 C = 1.5 for k — 0.75. There are four different 
motions: a circular-drift motion, a zigzag-1 motion, a 
zigzag-2 motion, and a straight-falling motion as we 
shall explain below. 

First of all, we consider the region 7 < 0, where a 
particle is motionless when the external force is ab- 
sent. When the external force g = (0, —5) is added 
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Figure 2: (Colour on-line) Trajectory of (a) the 
straight-falling motion, (b) the zigzag-1 motion, (c) 
the zigzag-2 motion and (d) the circular-drift motion. 
The arrows indicate the direction of migration. The 
parameters are chosen as 7 — 2 and K — 0.75 for (a) 
and 7 = 3 and k = 0.2 for (b), (c) and (d). 



and the magnitude is small, it causes a straight mo- 
tion in the direction of the external force as expected. 
Since b is chosen to be negative, the elongation is per- 
pendicular to the external force. We call this trivial 
motion a straight-falling motion. It should be em- 
phasized, however, that a non-trivial behavior oc- 
curs for large magnitudes of g as can be seen around 
g « 2 and 7 s» —0.01 in Fig. [TJa) where the straight 
motion is unstable and the zigzag-1 motion appears. 
The straight-falling motion and the zigzag-1 motion 
are displayed in Fig. [2ja) and Fig. [2jb), respectively. 
The bifurcation between the straight-falling motion 
and the zigzag-1 motion exists also for < 7 < j c 
as shown in Fig. HJa). When the external force is 
absent, this is the region that the particle under- 
goes a straight self-propelled motion, whose direction 
depends on the initial conditions. Figure |Tfb) for 
K = 0.75 indicates that such a bifurcation does not 
exist for 7 < 7 C . These results can be understood 
theoretically as described in section O 

Next, we consider the case 7 > 7c where a circular 
motion appears when g — 0. In the region indicated 
by the cross in Fig. []] where the magnitude of the 
external force g is Unite but small, a particle takes 
a circle-like motion. However, in this case, the cen- 
ter of the circle drifts to some direction as shown in 



Figure 3: (Colour on-line) Attractors in the ip-4> plane 
for (a) the straight-falling motion, (b) the zigzag-1 
motion, (c) the zigzag-2 motion and (d) the circular- 
drift motion. These are obtained by solving Eqs. (gj) 
- (JT]) numerically. The parameters for each motion 
are the same as those in Fig. [2] The arrows indicate 
the direction of motion. The square in (a) and the 
circle in (b), (c) and (d) indicate the stable and the 
unstable fixed points, respectively, and x indicates 
the saddle point. 



Fig. [DJd). We call this motion a circular-drift mo- 
tion. The drift direction asymptotically in time is 
determined uniquely for given values of the parame- 
ters. But it is neither equal to the direction of the 
external force nor determined by the initial condition. 
We discuss this property in detail in sectional For 
larger values of the external force, a zigzag-1 motion 
appears. For k = 0.2 in Fig. QJa), there is a coexist- 
ing parameter region of the circular-drift motion and 
the zigzag-1 motion. It is noted that there is another 
zigzag-like motion, called zigzag-2 motion, in the co- 
existing region, which is shown in Fig. HJc). In the 
case k = 0.75 in Fig. Hfb), the zigzag-2 motion has 
been observed only in the small region near g = 0.5 
and 7 = 4. When the external force is further in- 
creased, the zigzag-1 motion becomes unstable and a 
straight-falling motion appears for n — 0.75. This is 
contrast to the case of K = 0.2 where the straight- 
falling motion does not appear for 7 > j c . 

It is convenient to represent the motions asymp- 
totically in time in the (jj-ip plane where <f> and tp 
have been defined by (|8|) and ((9]), and (fT2|) respec- 
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Figure 4: (Colour on-line) The normalized frequency 
Ld/ujQ of the circular-drift motion (CD), the zigzag- 
1 motion (Zl), and the zigzag-2 motion (Z2) as a 
function of g for 7 = 3 and k = 0.2. 



tively. It is noted that the equations (J4J) - (jTj) are 
invariant under the transformations <fi — > <f> + 27r 
and 1^ — ^ ^ + 7T- Therefore, we may consider only 
the restricted range < <j) < 2tt and < V < 71 ■ 
The straight-falling motion is represented by the sta- 
ble fixed point (4>,ip) — (37r/2, tt/2), as shown in 
Fig. [31(a). The zigzag-1 and -2 motions are limit- 
cycles around this fixed point as shown in Fig. OUb) 
and (c) respectively. We discuss the difference of 
these motions in detail in the section[6] The attractor 
of the circular-drift motion is shown in Fig. |3jd). 

Now we discuss the coexistence region of the 
circular-drift motion, the zigzag-1 motion and the 
zigzag-2 motion in Fig. [IJa). Since the trajectories 
of these motions in Fig. [3] are periodic in the cfhip 
plane, one can define the period T for each motion. 
In Fig. 01 we display the frequency of these three mo- 
tions in the coexisting region for 7 = 3 and k = 0.2, 
which is normalized by the frequency luq of the cir- 
cular motion without the external forcing. We find 
that the frequency of the zigzag-1 motion becomes 
close to that of the circular-drift motion in the vicin- 
ity where the zigzag-1 motion becomes unstable. The 
frequency of the zigzag-2 motion is approximately a 
half of that of the zigzag-1 motion. 



4 Circular-Drift Motion 

One of the most interesting properties of the circular- 
drift motion is, as shown in Fig. [21(d), that the direc- 
tion of the drift motion is almost perpendicular to the 
external force g = (0, —g) with g > 0. The drift di- 
rection is defined quantitatively as follows. It should 
be noted that although the trajectory of a circular- 
drift motion in the real space is not periodic, it is 
precisely periodic in the velocity space. Therefore, 
we may define the position after one period starting 
from the coordinate origin as 



X 
Y 



dt v(t) 



d<f> 



COS <j){t) 

sin (p(t) 
-1 

«(0) 



COS ( 

sin < 



(16) 



where T is the period. The angle of the mean dis- 
placement with respect to the x-axis is given through 
the relation 

Y 

tan?7=^. (17) 

In Fig. [5j we show the angle 77 as a function of g 
for 7 = 3 and k — 0.75. Since the system pos- 
sesses the right-left symmetry, one may restrict to 
— 7r/2 < 77 < tt/2 without loss of generality. The 
angle rj is slightly smaller than 0, i.e., 77 f=s — 0.057T 
and gradually increases by increasing g up to the bi- 
furcation threshold g « 0.215. One unexpected phe- 
nomenon is that the value of rj becomes positive in 
the very vicinity of the threshold although the magni- 
tude is extremely small. This means that the particle 
gradually moves upward on an average while drifting 
to the right. Since the particle migrates by consum- 
ing the internal energy, this does not violate energy 
conservation. However, we do not have any definite 
explanation of this phenomenon which occurs in an 
extremely restricted parameter region. 



5 Numerical Results II 

In this section, we show the numerical results of self- 
propulsion in the electric-like external field. We have 
solved Eqs. (0J - © with g a = and Q afj ^ by 
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Figure 5: (Colour on-line) Angle of the mean dis- 
placement of the circular-drift motion as a function 
of the external force g for 7 = 3 and k = 0.75. 



using the fourth Runge-Kutta method with the time 
increment St ~ 1CP 4 for k — 0.5, a = —1.0 and 
b = -0.5. Therefore B = ab/2n = 0.5 > 0. 

Figure [5] displays, on the 7-/1 space, the phase di- 
agram of a variety of dynamical states: a circular 
motion, a zigzag-1 motion, a zigzag-2 motion, and a 
straight motion. Figure [7] shows the trajectory of the 
center of mass of each motion in the real space. We 
also show the trajectories of these motions in the 4>-ip 
plane in Fig. [8] When the magnitude of the exter- 
nal force h is large enough, a particle undergoes a 
straight motion for all 7. This occurs in the region 
indicated by the up and down triangles in Fig. [5] 
In this straight motion, even though we have chosen 
b < 0, the direction of the deformation is parallel 
to the external force in the region of the down tri- 
angles. The elongation becomes perpendicular for 
smaller values of h as indicated by the top trian- 
gles in Fig. [6] When the magnitude of the force h 
is decreased, the motion of a particle is affected by 
the three stable states which occur when the external 
force is absent [IS]: the motionless state for 7 < 0, 
the straight motion for < 7 < j c , and the circu- 
lar motion for 7 > j c . For 7 < 0, a particle does 
not move under a finite but weak external force, and 
therefore, the stable state is a motionless state. By 
increasing the magnitude h, a bifurcation occurs from 
the motionless state to a straight motion as shown 



Figure 6: (Colour on-line) Phase diagram on the 7-/1 
plane for k = 0.5 obtained by numerical simulations 
of Eqs. f4| - ([7|). The symbols indicate the circular 
motion (cross), the zigzag-1 motion (diamond), the 
zigzag-2 motion (square), and the motionless state 
(circle). The up and down triangles indicate the 
straight motion with n perpendicular and parallel 
to E respectively. The transition boundary between 
them is consistent with the analytical results given 
by Eq. which is shown by the thick dotted line. 
The thick solid line is the Hopf bifurcation bound- 
ary given by h = hn whereas the thin solid line is 
the pitchfork bifurcation boundary given by Eq. (|6T[) . 
The thin dotted line is the bifurcation boundary be- 
tween the circular motion and the zigzag-1 motion 
obtained numerically from the reduced equations (1471) 
and 



in Fig. [6] For < 7 < j c , a particle undergoes a 
straight motion irrespective of the values of the ex- 
ternal force. When 7 > j c , a particle undergoes an 
circular motion along an elliptically-deformed trajec- 
tory as in Fig. [T^d) under a finite but weak external 
force in the region indicated by the crosses in Fig. [6] 
Between this circular motion and the straight mo- 
tion, there is a region indicated by the diamonds in 
Fig. [5] where a zigzag-1 motion occurs as displayed in 
Fig. E^b). When 7 — j c > is large enough, there 
appears another motion, which is called a zigzag-2 
motion between the circular motion and the zigzag- 
1 motion as indicated by the square in Fig. [5] The 
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Figure 7: (Colour on-line) Trajectories of the center 
of mass in the real space for 7 = 3 (a) a straight 
motion for h = 1.3, (b) a zigzag-1 motion for h = 0.6, 
(c) a zigzag-2 motion for h = 0.5, and (d) a counter- 
clockwise circular motion for h — 0.3 obtained by 
solving Eqs. ((5) - © numerically. 
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Figure 8: (Colour on-line) Trajectories on the ip-cf> 
plane for (a) the straight motion, (b) the zigzag-1 
motion, (c) the zigzag-2 motion, and (d) the circular 
motion, obtained by solving Eqs. - CZJ) numeri- 
cally. The white circles and the white squares in- 
dicate the unstable and the stable fixed points, and 
the crosses indicate a saddle point. The values of the 
parameters are the same as those of Fig. [7] 



trajectory is displayed in Fig. [7[c) . 

It is evident in Fig. [6] that there is a region where 
the zigzag-2 motion and the zigzag-1 motion coexist. 
Although not clearly shown, there is a small region 
0.52 < h < 0.54 for 7 = 4 where the circular motion 
and the zigzag-2 motion coexist. These are also seen 
by analyzing the frequency of the motions. That is, 
the attractors of the circular motion and the zigzag-1 
motion and the zigzag-2 motion are periodic in the 
<j)-ijj space as shown in Fig. HI Therefore, one can 
define the frequency for these three motions, which is 
depicted in Fig. [9] where the frequency is normalized 
by that of the no external force ujq for (a) 7 = 1.5 
and (b) 7 = 4. Figure [^b) clearly indicates that 
there are coexistence regions of the circular motion 
and the zigzag-2 motion, and the zigzag-2 motion and 
the zigzag-1 motion. 

6 Analysis of the bifurcations I 

The set of time-evolution equations ((4]) - (J7J with 
Qa/3 — is complicated for theoretical analysis. Here 



we make a reduction of the variables by putting 
dv/dt = ds/dt = 0. This approximation is justified 
if the variables s and v relax rapidly compared with 
(j> and ip. Unfortunately, however, this is not gen- 
erally guaranteed in the present problem. Therefore, 
we need to check the results of the reduced system by 
comparing with the numerical results of the original 
system governed by Eqs. (0} - ©. 

Eliminating s and v from Eqs. (J4jl - dTj) , one obtains 

^=F(^) (18) 

^=G(0,VO, (19) 

where 

F ((f>, ip) = -Bv 2 cos 2^ sin 2ip - - cos (f> (20) 

v 

G(0,V) = ~ tan 2?/;- F (</>». (21) 

and B has been defined by Eq. (Tl4l . The magnitude 
s is given from Eq. ([6]) by 

s((t>,ip) = -v 2 cos 2ip. (22) 
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Figure 9: (Colour on-line) Normalized frequency of 
the circular motion (C), the zigzag-1 motion (Zl) and 
the zigzag-2 motion (Z2) for (a) 7 = 1.5 and (b) 
7 = 4. 



In the same way, the magnitude of the velocity for 
small value of g is determined from Eq. (|4]), as the 
largest positive root of the following cubic equation, 



3pv + 2g = 0, 



(23) 



where 



P(<t>,tp) = 



7 



3(1 + 5 cos 2 2ip) 

gsnicf) 
2(1 + 5 cos 2 2ip)' 



(24) 
(25) 



We have verified, by numerical simulations of the re- 
duced equations (fT8)l and (fl9|) . that all types of the 
motions except for the zigzag-2 motion are repro- 
duced for finite values of g. The dynamical phase 
diagram is quite similar to Fig. [T] obtained from the 
original equations (JU) - (0. 

If the variable s is also retained as a slow variable, 



one obtains 

lit 
ds 

~dt 
# 
~dt 



1 ■ 1 / 9 
-s sin zip cos 



-ks H 
bv 2 



bv cos 2ip 

.2 



as 



2s 



sin 2ip - 
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- COS( 
V 



(26) 
(27) 
(28) 



The magnitude of the velocity is given by the largest 
positive root of the cubic equation 



where 



p((p, ip) 



ipv + 2q = 0, (29) 

7/3- (a/6)scos2V> (30) 
(9/2) sin 4,. (31) 



It is verified numerically that these three-variable 
equations exhibit the zigzag-2 motion as well as other 
three kinds of motions mentioned above. However, 
since the analysis of the three- variable system is com- 
plicated, we carry out the stability analysis for the 
two-variable equations (fTS)) and (fH?j) below. 

Equations (fig)) and (fT^]) have two fixed points, 
(0,VO = (7r/2,7r/2) and (3tt/2, vr/2). We define the 
stability matrix as 



(32) 



where d^F — dF/dcf>. The stability matrix at the 
fixed point (7r/2,7r/2) is given by 



/7T 7T 

V2' 2 



9/v 

-g/v 



- 2Bv 2 
k + 2Bv 2 



(33) 



where v is a positive root of Eq. (|23[) with p — pq and 
q = qo where 



<7o = 



7 



3(1 + 5) 
g 

2(1 + 5) 



(34) 
(35) 



The trace and the determinant of the stability matrix 
are calculated as 

tri(vr/2,7r/2) =g/v - k + 2Bv 2 (36) 
det L (tt/2, tt/2) = -gn/v. (37) 
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Since v should be positive, the determinant is always 
negative so that we find the fixed point (it/2, tt/2) is 
a saddle point. 

The stability matrix for the other fixed point 
(37r/2,7r/2) corresponding to the straight-falling mo- 
tion is given by 



3tt 



2 ' 2 

from which one obtains 



-g/v 
g/v 



- 2Bv 2 
-K + 2Bv 2 



(38) 



(39) 
(40) 



tr .L(37r/2,7r/2) = -g/v- n + 2Bv 2 
detL (37r/2,7r/2) =gn/v. 



where v is a positive root of the cubic equation (1231 
with p = Po and q = —qa- Since the determinant 
(j4H)l is always positive, a Hopf bifurcation occurs at 
trL(37r/2, tt/2) = 0. At this fixed point, Eq. can 
be written as 



(1 + B)v 3 - jv 



0. 



(41) 



Using this equation, the derivative of the trace (|39|) 
with respect to g is calculated as 

d T /3tt tt\ 2(B- 1)v 2 
— trL — - = —i —L . 42) 

dg \ 2 ' 2J 2(1 + B)v 3 +g V ' 

Therefore, the stability condition of the fixed point is 

(43) 



>g' 
<g* 



for B < 1 
for B > 1 



where the bifurcation threshold g* is determined as 
follows. From Eqs. ^ and dUJ, v(> 0) at the 
threshold is obtained as 



7 — K 
1-B 



1/2 



(44) 



provided that (7 — k)/(1 — B) > 0. Substituting v 
into Eq. (|41[) . one obtains the bifurcation boundary 
as 
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2£( 7 - 7c ) ( 1 



1 - B 



1/2 



(45) 



(1-B) 

Since g* > 0, the following condition should be sat- 
isfied; 

B < 1 for 7 > 7c 
£? > 1 for 7 < 7 C 



(46) 



The bifurcation threshold (l4"5j) is indicated by the 
thick solid line in Figs. [TJa) and (b). Note that 
B = 1.25 for k = 0.2 and B — 1/3 for = 0.75. 
The stability condition (|43"1) and the bifurcation lines 
agree with the numerical results obtained from the 
original equations (JU) - (0. 

When 7 > 7c, there exists another bifurcation 
from the circular-drift motion to a zigzag- 1 motion 
as shown in Figs. HJa) and (b). We have obtained 
the boundary of this bifurcation by numerical simu- 
lations of the reduced equations (fT8|) and (IT9l) . The 
results are plotted by the thin solid lines in Figs. [II a) 
and (b). The line for n = 0.75 is in an apparent 
agreement with the numerical results of the origi- 
nal equations (@} - ©• See, however, further dis- 
cussion given in the next paragraph. The thin solid 
line for k = 0.2 in Fig. [Ua) a l so ag ree s with the sta- 
bility threshold of the circular-drift motion but the 
reduced equations do not reproduce the coexistence 
of the zigzag-1, zigzag-2 and circular-drift motions. 

We have examined the motion for k = 0.75 in the 
vicinity of the bifurcation between the circular-drift 
motion and the zigzag-1 motion. The trajectory in 
the real space obtained numerically from the original 
equations (0J - © for 7 = 3 and g = 0.232 is dis- 
played in Fig. [Tola). It is clear that this is neither 
a simple circular-drift motion nor a zigzag-1 motion 
but is a kind of a mixture of these two motions, or a 
mixture of a circular-drift motion and a zigzag-2 mo- 
tion. In fact, the trajectory in the tfi-tp space shown 
in Fig. [TOTb) is found to be a superposition of the 
circular-drift motion and the zigzag-2 motion. Sec 
Figs. Etc) and (d). The bifurcation for the reduced 
equations (|T8l) and ([19]) is at about g = 0.245 for 
k = 0.75 and 7 = 3. The trajectories near this thresh- 
old are shown in Fig. [10]for (c) g = 0.245132 and (d) 
g = 0.245133. There are two trajectories in Fig. fTOt c') 
corresponding to a right-moving and a left-moving 
circular-drift motions but these two merge each other 
near the saddle point in Fig.fTUTd) to cause a zigzag-1 
motion. Therefore this is a saddle homoclinic orbit 
bifurcation. The bifurcation exhibited in the original 
set of equations is more complicated. 
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(a) 



(c) 










o 





% 2% % 2k 

Figure 10: (Colour on-line) Trajectory in the vicinity 
of the bifurcation between the circular-drift motion 
and the zigzag-1 motion for 7 = 3 and k = 0.75 (a) in 
the real space and (b)-(d) in the 0-0 space, (a) and 
(b) are obtained from the original equations (J4]) - (JTJ) 
for g = 0.232 whereas (c) and (d) are obtained from 
Eqs. HHJ and (31) for g = 0.245132 and 0.245133 
respectively. The cross and circle in Fig. (b) - (d) 
indicate the saddle point and the unstable fixed point, 
respectively. 



7 Analysis of the bifurcations 
II 

In this section, we analyze the bifurcations obtained 
by the original equations (J4j) - ((Tj) with 5 = 0. As in 
the preceding section, we consider a simplified set of 
equations eliminating s and v. From Eqs. Q - ([7]), 
by putting dv/dt = ds/dt = 0, we obtain 



with 



G(0,0) 



F((f>,tp) = — sin20 
bv 2 - 



(47) 
(48) 



(49) 



2s 



■ sin 2-0 



2s 



sin 2(0 + ?/>), 



where s is given by 

s(0, 0) = —v 2 cos 2tp + — cos 2(0 + ijj) 



(51) 



and v is the largest positive root of a cubic equation 
(1 + Bcos 2 2ip)v 3 -Tv = 0. (52) 
We have defined 

r(0, tp) = 7 - (h/b)B cos 2(0 + 0) cos 20. (53) 
From Eqs. (j52")l and (|5"Tj) , v and s are given for T > 

by 



l! ((/), 1p) 
S (0, 1p) = 

and for T < by 



r(0,VO 



1/2 



1 + i? cos 2 2^, 
67 cos 2ip + h cos 2(0 + -0) 
k(1 + BCOS 2 20) 



«(0,0) =0 
,0) = (/i//«)cos2(0 + V). 



(54) 
(55) 



(56) 
(57) 



In both cases, s > is required. By the numerical 
simulations of the reduced equations (T4"T|) and (|4"5|) . 
all types of motions except for the zigzag-2 motion 
are obtained. 

As in section [BJ we have verified numerically that 
the zigzag-2 motion is realized if we retain s as a slow 
variable; 



ds 

— = — KS 
dt 

dip bv 2 
~dt ~ 



+ &v 2 cos20 + /icos2(0 + 0) (59) 
— (xs^ h 

sin 20 - — sin 2(0 + ip) (60) 



where 



(50) 



f (7 - §scos2?/>) 1/2 if 7- fscos2V> > 
V ~\ if 7- §scos2V> < 

(61) 

However, in the theoretical analysis given below, we 
employ Eqs. (|47|) and (l48l) which are much easier to 
treat. 
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In the restricted space < (f) < 2ir and < ip < ir, 
Eqs. (|47l) and (|48l) have four fixed points, (0, ip) — 
(n7r/2,^ ) with n = 0, 1,2,3. In order to satisfy the 
condition s > 0, we note from Eqs. and ([57)1 that 
rpo = tt/2 for n = 1 and 3 while for n — and 2, -0o = 
when 67 + ft, > and V^o = tt/2 when 67 + ft < 0. 
If r = &7 — ft£> > 0, i.e. v > 0, the two solutions 
(<fi,ip) = (mr/2, ipo) with 71 = and 2 represent a 
straight motion propagating to the right and the left, 
respectively, whereas, if T < 0, i.e. v — 0, they are a 
pair of degenerate solutions of the motionless state. 

We introduce the linear stability matrix in the 
same form as Eq. (|32|) and analyze the stability of 
the four fixed points. In the region r > 0, v and s 
are given by Eqs. ([M]) and (jB1)j) respectively. From 
the stability matrix around the fixed points (O,'0o) 
and (71", i/jo) j the trace and the determinant of L((j), ip) 
are given, respectively, by 

2B(7 - 7 C ) + ah/n 



tvLi + = 



det L\ + = —ah. 



(62) 
(63) 



When a < and ft, > as we have assumed, the trace 
is positive if < ft < hn = |fo|(7 — 7c) for 7 > 7 C and 
negative if ft, > hu for 7 > 7 C or for 7 < 7 C , whereas 
the determinant is always positive. Therefore, the 
fixed points (0,tpo) and (tt,^q) corresponding to a 
straight motion are always stable if 7 < 7 C as ex- 
pected. When 7 > 7c, they are stable for ft > hn and 
they lose their stability for < h < ft# by a Hopf 
bifurcation. This corresponds to a transition from a 
straight motion for h > hn to a zigzag- 1 motion for 
< h < h H - See Figs. UJa) and (b). The bifurca- 
tion boundary hn represented by the thick solid line 
in Fig. [6j is found to agree with the results of the 
numerical simulations of the original equations ((4]) - 
©• 

There is another bifurcation in the region where 
the fixed points (O,-0o) and (tTjV'o) for the straight 
motion are stable. The requirement s > leads from 
Eq. ([55]) to 



tt/2 for ft < ft* = |6|7 

for ft > ft* 



(64) 



This means that elongation of the particle is perpen- 
dicular (parallel) to the velocity (external force) for 




Figure 11: (Colour on-line) Migration velocity in 
the vicinity of the pitchfork bifurcation boundary 
7 = —0.5. The dots indicate the results obtained 
by solving Eqs. (HJ) - ((TJ) numerically and the dotted 
line indicates the velocity obtained analytically from 
Eq. (JSHJ). 



ft < (>)ft*. Since we have assumed ft > 0, this oc- 
curs only for 7 > 0. This bifurcation threshold repre- 
sented by the thick dotted line in Fig. [6] is consistent 
with the results of numerical simulations of the orig- 
inal equations (|4j) - ([TJ) . 

The determinant and the trace of the stability 
matrix of the other two fixed points (7r/2,-0o) and 
(37r/2,?/>o) are given, respectively, by 



trio 



= 2B(j - 7 c) - ah/n 

(1 + B) 
det L2+ — ah. 



(65) 
(66) 



The determinant is always negative for a < and 
ft > 0. This implies that the fixed points (7r/2, V»o) 
and (37t/2,-0 o ) are saddle points. 

In the region T < 0, we can show by a similar 
analysis that the fixed-points (0, ^0) and (77, ^o) with 
'ipo = are stable, whereas the other two fixed-points 
(n/2,'tp ) and (37r/2, i/i ) are saddle. 

The bifurcation from the motionless state to the 
straight motion in Fig. [6] is derived as follows. For 
the stable straight motion (0, VO = (0,0) and (7r,0), 
T is given by V — 7 - (h/b)B. We note from Eq. ([52]) 
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that a bifurcation from the state v — to v ^ 
occurs at 7 = (h/b)B. Since h > 0, B > and 6 < 0, 
this bifurcation is possible only for 7 < 0. Therefore, 
when h is increased, the bifurcation occurs at h = h p , 
where 

15,7 (67) 



h p = 



B 



This bifurcation boundary indicated by the thin solid 
line in Fig. [5] agrees quite well with the results of 
numerical simulations of the original equations ((4]) - 
([7|). The velocity around the pitchfork bifurcation is 
given from Eq. ([54)) by 



v pf 



j + hB/\b\ 
l + B 



1/2 



(68) 



This theoretical result is plotted in Fig. [TT] in com- 
parison with the numerical result of Eqs. (U) - (0. 

Finally, we mention that the thin dotted line in 
Fig. [6] is the boundary between the circular motion 
and the zigzag- 1 motion obtained numerically from 
Eqs. (|47|) and (|48|) . Since this set of reduced equa- 
tions does not reproduce the zigzag-2 motion, the 
agreement with the results of the original equations 
(gl) - is poor. 

To summarize, the reduced equations (|47|) and (|48| 
provide very accurately the three bifurcation thresh- 
olds given by the thick solid line, the thick dotted line 
and the thin solid line in Fig- El However, in order to 
reproduce theoretically the zigzag-2 motion, we have 
to take into consideration of the variable s. 

8 Discussion 

We have investigated dynamics of a self-propelled 
particle under two types of external forcing, the 
gravitational-like force and the electric-like force. In 
both cases, we have found a variety of dynamical mo- 
tions and have obtained numerically the dynamical 
phase diagrams (Figs. [T]and|n]). We have also ana- 
lyzed the bifurcations of these motions by using the 
reduced equations in terms of the two kinds of an- 
gles which represent the migration velocity and the 
elongation of a particle. 

In the case of the gravitational-like force, we have 
found the circular-drift motion, the zigzag- 1 motion, 



the zigzag-2 motion, and the straight-falling motion. 
The reduced equations (fT8f and (H"9|) admit these 
solutions except for the zigzag-2 motion. We have 
shown that a Hopf bifurcation appears between the 
zigzag- 1 motion and the straight-falling motion. This 
property agrees with the numerical results of the orig- 
inal equations ([TJ and ([2]) with Q a fj = as shown in 

Fig. m 

In the case of the electric-like force, we have ob- 
tained the circular motion, the zigzag- 1 motion, the 
zigzag-2 motion, and two types of the straight mo- 
tion, whose direction of deformation is either parallel 
or perpendicular to the velocity vector. The reduced 
equations (|4"T|) and (|4"5f reproduce all of these motions 
except for the zigzag-2 motion. It has been shown 
that there appears a Hopf bifurcation between the 
zigzag- 1 motion and the straight motion, and a pitch- 
fork bifurcation between the motionless state and the 
straight motion. These properties are consistent with 
the numerical results of the original equations (fTJ) and 
([2]) with g = as shown in Fig. [6j In the both cases, 
however, the variable s has to be considered as a slow 
variable to realize theoretically the zigzag-2 motion. 

The basic ingredient of the systems studied is the 
competition between the external forcing and the cir- 
cular motion which occurs in the absence of the ex- 
ternal forces. This causes the transitions between 
the straight motion and the two types of zigzag mo- 
tions and the dynamics of circular-drift motion. At 
present, we have no experiments which are directly 
related to the predictions. However, since our time- 
evolution equations are general based on the symme- 
try argument, we expect that the predictions made 
in the present paper will be detected experimentally 
in the near future. 
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